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Abstract 

A radial basis function implementation of the 
meshless local Petrov-Galerkin (MLPG) method is 
presented to study Euler-Bernoulli beam problems. 
Radial basis functions, rather than generalized moving 
least squares (GMLS) interpolations, are used to 
develop the trial functions. This choice yields a 
computationally simpler method as fewer matrix 
inversions and multiplications are required than when 
GMLS interpolations are used. Test functions are 
chosen as simple weight functions as in the 
conventional MLPG method. Compactly and non- 
compactly supported radial basis functions are 
considered. The non-compactly supported cubic radial 
basis function is found to perform very well. Results 
obtained from the radial basis MLPG method are 
comparable to those obtained using the conventional 
MLPG method for mixed boundary value problems and 
problems with discontinuous loading conditions. 

Introduction 

Meshless methods are developed to overcome 
some of the disadvantages of the finite element method 
(FEM) such as discontinuous secondary variables 
across inter-element boundaries and the need for 
remeshing in large deformation problems. 1 " 4 Recent 
literature shows extensive research on meshless 
methods and, in particular, the meshless local Petrov- 
Galerkin (MLPG) method. The majority of literature 
published to date on the MLPG method presents 
variations of the method for C° problems. 5 ' 6 However, 
a comparatively limited amount of work 4 7-10 is 
reported on the more complicated C 1 problems. Atluri 
et al. 4 present an analysis of thin beam problems using 
a Galerkin implementation of the MLPG method. In 
reference 4, a generalized moving least squares 
(GMLS) approximation is used to construct the trial 
functions, and the test functions are chosen from the 
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same space. In references 11-14, a meshless Petrov- 
Galerkin implementation of the MLPG method is 
presented; the GMLS approximation is used to 
construct the trial functions, and the test functions are 
chosen from a different space. Closer scrutiny of these 
formulations shows that a large number of calculations 
are required to compute the first and second order 
derivatives of the moving least squares (MLS) trial 
functions. Hence, a computationally efficient 
alternative to the MLS trial functions is preferred. 

This paper demonstrates the use of radial basis 
interpolation functions in the meshless local Petrov- 
Galerkin formulation for beam problems. The radial 
basis functions are simple, and the evaluation of the 
derivatives is simpler than for the traditional MLS 
approximations. In the present radial basis MLPG 
(RPG) formulation, simple weight functions are chosen 
as test functions, and Gaussian quadrature is used to 
integrate the weak form. The effectiveness of the RPG 
method is evaluated by applying the formulation to a 
variety of patch test and mixed boundary value 
problems. 

The outline of the paper is as follows: First, the 
moving least squares interpolation used in the 
conventional MLPG method is discussed as motivation 
for finding a more computationally efficient alternative. 
Next, an overview of radial basis functions (RBF) for 
C° problems is presented; the shape functions obtained 
from radial basis interpolation are derived, and the 
shape functions obtained when polynomial basis 
functions are included in the interpolation are derived. 
The development of these radial basis shape functions 
is then expanded and repeated for beam problems. The 
system of algebraic equations developed from the local 
weak form of the governing differential equation and 
the chosen trial and test functions is presented. Patch 
test problems are used to validate the RPG method for 
different choices of radial basis function. Then, the 
RPG method is applied to mixed boundary value 
problems. Finally, the method is applied to a problem 
with discontinuous loading conditions. 

Interpolation Schemes 

In this section, the moving least squares 
interpolation scheme used in the conventional MLPG 
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method is discussed first. Then, two interpolation 
schemes involving radial basis functions (RBF) are 
presented. In the first scheme, radial basis functions 
alone are used to construct the shape functions. The 
second scheme is a hybrid that uses both radial basis 
functions and polynomial basis functions to construct 
the shape functions. 

The Moving Least Squares Interpolation 

A moving least squares (MLS) interpolation is a 
scheme that passes a smooth function through an 
assumed set of fictitious nodal values. The 
interpolation is performed such that the least squares 
error between the function and the nodal values is a 
minimum. 1 ' 2 A schematic of the MLS interpolation is 
presented in Figure 1. 
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and A/x), j = 1 .../ 7 , is a weight function. The first 
derivative of these shape functions is all that is required 
by the MLPG method and is given by 13 
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C 1 MLS Shape Functions : The C 1 GMLS shape 
functions for deflection, w, and slope, 6, in 1-D are 
given, using the local coordinate approach of references 
12 and 14, by 


Figure 1: Moving Least Squares (MLS) interpolation 

C° MLS Shape Functions : In one-dimensional 
problems, the C° MLS shape functions are given by 
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In Equation (2), [P] is an («, m) matrix, and [X] is a 
diagonal (/?, n) matrix defined as 
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In Equations (8) - (10), p(x) is a polynomial basis 
function, and 
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In Equation (11), [P] and [PJ are (n, m) matrices, and 
[A,] is a diagonal («, n ) matrix defined as 
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where = x^ - x j , k= 1,2, ...,/?, 

p T (|) = [l, | 2 , ... J, and (14a) 


are given in reference 13 and are not repeated here. 
Note how the additional degree of freedom (9) and the 
need for the higher order derivatives yield very 
complicated expressions for these derivatives. For thin 
plate problems (2-D C 1 problems), these derivatives 
become even more complicated. (Expressions for the 
partial derivatives for the 2-D shape function may be 
found in reference 10.) Therefore, a more 
computationally efficient method for approximating the 
trial functions in the MLPG method is sought. Radial 
basis functions appear to be a good candidate for 
achieving such a purpose because the shape functions 
obtained from radial basis interpolation are simpler than 
the shape functions presented above for the MLS. 
More importantly, the derivatives of the radial basis 
shape functions are simple and involve considerably 
fewer matrix inverse and multiplication operations than 
the derivatives of the MLS shape functions. The radial 
basis functions (RBF) are discussed next. 
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In C 1 problems, the first, second, and third derivatives 
are required by the MLPG method. The first 
derivatives of y/j are 




(w) 

= n * = 


m ( / 

[BJ 

g = 1 


'gl 


+ 



+ [A]“ 1 [B W ] 



(16a) 


and 

= zk«([ A r‘[ b»]L- 

g=l 

+ ^([Ar 1 [B 0 ] >JC +[A]- 1 [B <? ]y 


(16b) 


where 


Radial Basis Function 


The radial basis formulation provides a continuous 
interpolating function for u(x) as a linear combination 
of radial functions. 15 The interpolating function is 
given by 

N 

u(x)= Y J R j(x)a j , (17) 

7=1 

where R/x). the radial basis functions (RBF), are 
functions at each of the N scattered points, and cij are 
the unknown coefficients, j = 1, 2, ..., N. The RBF, 
R/(x). are functions of distance r j and are defined as 

R j (x) = R j (r j ). (18) 

The radial distance, rj , in Cartesian coordinates can be 
expressed as 

r j(X) = xl{ x - x jf + {y-yjY , (19) 

where Xj and y, are the coordinates of node j. Equation 
(17) can be written in matrix form for C° problems as 

w(x) = R T (x)a , (20) 


[A]" 1 =4A]" 1 [A] >JC [A]- 1 . (7) 

The second derivatives and third derivatives involve 
considerably more complex expressions containing 

[A]” 1 , [A]" 1 , [A]^, etc., and the detailed derivations 


where 
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Forcing the interpolation of Equation (17) to pass 
through the N scattered points, a set of equations to 
determine the coefficients a, can be written as 


N 
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where 


(22) at the nodes exactly. As the number of non-nodal 
interpolation points, M, is increased, the shape 
functions in Equation (28) satisfy 
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Note that R B is an ( N , N) matrix. Here, 
u ( uj , j = 1, 2, ..., N) are the nodal values of u at the 

N scattered points. The unknown coefficients in 
Equation (22) can be obtained as 


M 
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Unlike in the moving least squares (MLS) method, 
the derivatives of the shape functions are easy to 
evaluate using Equation (28) as 
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The interpolating function for m ( x ) in Equation (20) can 
now be rewritten as 
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where 4i/ are the elements of the matrix R B . The 

shape function <pj(x) obtained through the above 
procedure satisfies the Kronecker Delta property only at 
the nodes \ i.e., 

<Pj( x j) = 1 , and 

(29) 

<Pj( x k) = 0- 

Note that the shape functions in Equation (28) also 
satisfy the property 


Some of the classical radial functions used in 
multivariate interpolation are presented in Table 1. 
Note that the shape parameter c in the radial basis 
functions in Table 1 is user-defined and can be adjusted 


Table 1: Classical radial basis functions 16 


Classical RBF 

Equation 

Linear 

R(r) = cr 

Cubic 

R(r) = (/' + <?) 3 

Thin plate spline 

R(r) = r 2 log(c/' 2 ) 

Gaussian 

R(r) = e~ cr2 

Multiquadric 

R(r) = -Jr" +c“ 
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(39) 


to fit the required data. The classical radial functions 
have two limitations; (1) the matrix R B may not be 
positive definite, and (2) the functions do not possess 
local support, i.e, changing the location of the center 
( Xj,y>j ) in Equation (19) affects the entire interpolation, 
to overcome these limitations, compactly supported 
positive definite radial functions were proposed. 17 
These functions were derived using a constraint to 
guarantee positive definiteness of the interpolation 
matrix, R B . The compact support of these functions 
guarantees that every point in a compact radial basis 
interpolation domain does not necessarily have an 
affect on every other point in the domain. The compact 
RBF adapted and used here are 

Compact-I: 


f = l + x, 

respectively. The function values are evaluated at the 5 
nodes and are prescribed as u/s in Equation (26). The 
values of u are evaluated using Equation (26) at 200 
points in the interval x\ < x < x 2 and are plotted in 
Figures 2b and 3b. As seen from these figures, the 
compact RBF recovers the polynomial values 
(Equations (38) and (39)) exactly only at the 5 nodal 
points, and elsewhere the values of the polynomials in 
Equations (38) and (39) are not recovered. 

In order to improve the polynomial accuracy of the 
solutions, Powell 15 suggested adding polynomial basis 
functions to the radial basis functions as 
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where (t = r / sj), and Sj is the radius of the domain of 
compact support. The shape functions (Equation (28)) 
obtained from the Compact-I and Compact-II functions 
possess all the properties in Equations (29) - (31). 
Several other forms of compact support functions can 
be found in references 16 and 17. 

Hybrid Radial Basis Function 



The classical radial basis functions shown in Table 
1 and the compactly supported functions in Equations 
(36) and (37) cannot represent polynomial solutions 
exactly 18, 19 ; they can represent the polynomial values 
only at the N scattered points. Figures 2 and 3 show 
radial basis interpolations obtained from the compact 
RBF in Equation (36) with Sj = 0.6 using 5 nodes in the 
interval <x<x 2 , where X\ = -1 and Xi = +1. 
Figures 2a and 3a show the RBF values that correspond 
to a constant polynomial, 

/ = 1 > (38) 

and a linear polynomial. 


N m 
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where R h cij, and N are as in Equation (17), p(x ) is the 
polynomial basis function, z k are the unkown 
coefficients associated with the k th polynomial term, 
and m is the order of the polynomial basis function. 
Equation (40) is written in matrix form as 18 


u(x) = R T (x)a + p T (x)z 
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where a and R T (x) are as in Equation (21), and 
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The interpolation is forced to pass through the N points, 
with the constraint, 

N 

YjPk(Xj) a j =° , k=l,2,...,m, (43) 

./=! 

imposed to guarantee unique approximation. 18 The set 
of equations to determine the coefficients cij and z k is 
thus written as 
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Node: 1 2 3 4 5 

Figure 2a: Radial basis function values that correspond 
to a constant polynomial 



Node: 1 2 3 4 5 

Figure 2b: Interpolated and exact values of a constant 
polynomial 



Node: 1 2 3 4 5 

Figure 3a: Radial basis function values that correspond 
to a linear polynomial 



Node: 1 2 3 4 5 

Figure 3b: Interpolated and exact values of a linear 
polynomial 


where u r and R B are defined in Equations (23) and 
(24), and 
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The unknown coefficients in Equation (44) are obtained 
as 


The nodal shape functions are then 

N m 

<Pl (x) = Z R j (x)X jl + Z Pk ( x )r(N+k)(l) , (48) 

7=1 *=1 

where ;/ are the elements of the matrix [G]" 1 . The shape 
functions in Equation (48) possesss the Kronecker 
Delta and the unity partition properties of Equations 
(29) and (30). The derivatives of the shape functions 
are easy to evaluate using Equation (48) as 
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The interpolating functions m(x) in Equation (41) can 
now be rewritten as 
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where (dR j / dx) and (dRj / dy) are as in Equations 
(33) -(35). 
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Beam Problem Interpolation Schemes 


Similarly, Equation (52) is written as 


This section presents the interpolation schemes 
used in the RPG method for beam problems. The shape 
functions for both the radial basis and hybrid 
interpolations are derived. These shape functions will 
be used in the next section in the system of algebraic 
equations developed from the local weak form of the 
governing differential equation. 

Radial Basis Function 


The radial basis functions, R/(x), are functions of r h 
where in 1-D 

rj(x) = x-Xj . (50) 

In C 1 problems, the deflection, w, and the slope, 
9=dw/dx, are both primary variables and degrees of 
freedom whose continuity need to be satisfied. The 
interpolating function for w(x) is assumed to be of the 
form 


0 = 
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Forcing the interpolations to pass through N nodal 
values, the set of equations to estimate the coefficients 
cij and bj is written as 

[Q B ] W = {d} , (57) 

(2JV, 2JV) (2N, 1) (2N, 1) 

where 


w(v) = R\ (x)fli + 5) (x)b\ + R 2 ( x)a 2 + S 2 (x)b 2 
+ . . . + R N (x)a N + Sff {x)bft 
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where cij and bj,j= 1, 2, N, are unknown coefficients, 
Rj(x) are the radial basis functions, and Sj(x) = dRpc) / 
dx. Because of the direct relationship between the slope 
and the deflection, the approximating functions for 9 
cannot be chosen independently from the functions for 
w, and as in Equation (51), the approximations for #are 
written as 
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The unknown coefficients in Equation (57) are obtained 
as 


In matrix form, Equation (51) is written as 
w(x) = Q r (x){c} , 
where 

Q T (x) = [R l (x) S } (x) R 2 (x) S 2 (x) ... 

R n (x) 5 ^(x)] 
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The interpolation for w in Equation (53) can be written 
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where <p are the nodal shape functions, 


<p( x ) = Q t (^)[Q b ] _1 

= y[ w \x) y/[ 9 \x) y/ ( 2 V \x) yx {9 \x) ■■■ (62) 
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where {c} and Q T (x) are as in Equation (54), and p T (x) 
are the polynomial basis functions, 


P T (*) = [a(*)> P 2 (x), P3(x), ..., p m (x)\ 


II X x 2 ... x m ~ l 


(68) 


From Equation (62), the individual shape functions for 
deflection and slope, y/^’\x) and y/^ d \x),are, 
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where rjki are the elements of the matrix [Q B ] The 
derivatives of these shape functions are easily evaluated 
as 
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and {z} are the unknown coefficients associated with 

p t M, 

{z} = {zi, z 2 , z 3 , ..., z m j T . (69) 


The interpolation of Equation (67) is required to pass 
through the N points with constraints, 
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where k = 1, 2, ..., m, and cij and bj are the unknown 
coefficients in Equation (54). Equation (70) is imposed 
to guarantee a unique approximation. The set of 
equations to determine the coefficients {c} and {z} is 
thus written as 
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Hybrid Radial Basis Function 


As discussed for C° problems, in order to improve 
the polynomial accuracy of the solutions, interpolations 
involving both radial basis functions and polynomial 
basis functions are considered as 


w(x) = Q T (x) {c} + p T (x){z} 
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where {d} T and [Q B ] are defined in Equations (58) and 
(59), and [T B ] is a (27V, m) matrix, 
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The unknown coefficients in Equation (71) are obtained 
as 



(73) 


The interpolating function for vv in Equation (67) is 
now written as 


w(x) = 


|Q t « 


P T W 


I G\~ l 



2 N 

= Z <Pj( x )dj, 
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where <p are the nodal shape functions 


(74) 


subjected to four boundary conditions, two at each end 
(x = 0 and x = L). The boundary conditions are on vv, 
9 V, and M, where 

3 2 

_ dw rr „ T d vv , ,, „ T d vv 

9 = — , V = -EI — and M = El — — (78) 

d x dx 3 dx~ 

are the slope, shear force, and moment, respectively. 
The essential boundary conditions are on vv and 9, 
while the natural boundary conditions are on V and M. 
The boundary condition sets on vv and V and 9 and M 
are disjoint, i.e., if vv is prescribed then V cannot be 
prescribed, and vice versa. 

The MLPG equations are derived using a weighted 
residual weak form of the governing equation (Equation 
(77)). The MLPG equations are 4,11 ’ 13,14 


?>(x) = [q t M P T (x)J[GT 


= \//[ w) (x) y/[ 9 \x) y/ (9 \x) ... 

yX { N’\x) Yn\ X ) 


K (node) d + K (bdry) d _ f (node) _ f (bdry) _ q 

(75) where the superscript “bdry” denotes boundary, 

d T = {h’i 0\ w 2 0 2 ■■■ w N 0 N } (80a) 


From Equation (75), the shape functions for the 
deflection and slope, and \// {9 \x ) , may be 

written as 


are the nodal values of deflections, vv, and slopes, 9, at 
all the N nodes of the model used to analyze the 
problem (Equation (58)), and 


V'l " 1 (■ x) - Z { R j ( x )'£(2 /-1)(2/-1) + SjW'f(2 j)(2l-l) ) 

7=1 


+ Z Pk ( x ) ' f(2JV+i)(2/-l) 

k=\ 


- Z( S 7( x )'^"(27-l)(2/) + Sj( x )'C(2j)(2l)) 
7=1 


(76) 


+ Z Pk ( x ) ■ ({2N+k)(2l) 

7=1 


where £*/ are the elements of the matrix [ G ] _1 . The 
derivatives of these shape functions are easy to evaluate 
as in Equations (64) - (66). 

MLPG Equations for Beam Problems 

In this paper, the radial basis function is used in the 
MLPG method for beam problems that are governed by 
the fourth-order equation 

d 4 

EI—^ = f in 0 <x<L (77) 

dx 4 


K (node) =[ k | / node) J (80b) 

K (bdly) =[k| bdry) j (80c) 


with 




(node) _ 


= EI 


,2 (w) .J 2 (w) ,2 (w) J 2 „.(7b 

f dX \ d X f dX < d¥ i dx 

J (i \ dx dx dx dx 

dx- dx~ L dx~ dx 2 

ny' nyl 


+ n r EI 


dx~ dx 


Xi 


(iv) 


d y/ 




dx 5 


jV ( w) 

(0) u r ] 


dx 5 

m d M 6) 


Jr $ 


-n Y EI 


dx dx 

-7 ( w ) i (w) ^2 ...(0) 

dx } d Vj dx) d Vj 


dx dx 2 dx dx 2 

dtP dV? d X f> 1 


dx dx 2 


dx dx~ 


-(0 


(80d) 
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r (bdry) . 


x'rvp x'rvp 


x?\f\ T u 


+ n J'J 


„(»-) 


rfV <w) 

" y'j iw) a tj 


xi 


(0) 


dx i 

dVi w) 


dx 


Xi 


W) 


dx i 

d\f 


dx 


dO 



(80e) 


p(node) 


\x\ w) f 


dx 


\x\ 6) fdx 


and 


(80f) 



d X A 


r (») 

Ai 

f <bdry, =H r M< 

dx 

+ n x V- 


dxf ] 


y8» 


dx 

r (0 

Ai 


r (0 

sV 


V w) l 

Ai 


dxA 

• +a e 8< 

dx 


r (0 


{. Xi J 

1 SW 

dx 


AO 

S0 


(80g) 


where i = 1, 2, ..., N and j = 1, 2, n, and n is the 
number of nodes in the domain of definition of the trial 
function. In these equations, j/ w) and are 

components of the test functions, and y/j 0) are the 
shape functions, Q s (,) (see Figure 4b) is the local sub- 
domain of the test function at node i, n x is the unit 
outward normal to Q/‘ ) , and F s (,) are the boundary 
points of Q iS (,) (see Figure 4b). When Y® coincides 
with an interior point, that point is denoted T s i'\ and 
r sw (,) , r sM u \ and rj‘> denote the boundary points 
where F/' 1 intersects the boundary when w, 9, M, and V 
are prescribed, respectively. Also in these equations, 
a u and u 0 are penalty parameters to enforce the 
essential boundary conditions, and w, 8 , M , and V 
are prescribed values of the deflection, slope, moment. 


and shear, respectively. See reference 14 for a more 
detailed explanation of these terms. 

• • » t • m — *-x 

12 k j i N 

(a) An TV-node model of a beam 



(b) Components of the trial and test functions 


Figure 4: Comparison of the domains of the trial and 
test functions 


The system of equations presented in Equations 
(79) - (80g) is the general set of equations valid for any 
set of trial and test functions. In this paper, a Petrov- 
Galerkin method is used; the test functions are chosen 
to be different from the trial functions. The choices for 
the trial and test functions are now briefly discussed. 

Trial Functions 


The trial functions are chosen as 

Mx) = +y,f «*,•), (81) 

7=1 

where t^J w) (x) and are the radial basis shape 

functions of Equation (63), and Wj and 9 ; are the nodal 
values of w and <9 at the N nodes (Equation (58)). Note 
that in MLPG algorithms employing the moving least 
squares interpolation scheme for the trial functions, the 
values d (Equation (80a)) are fictitious nodal values, d . 
In this paper, radial basis functions are fit to the actual 
nodal values, d. 

Test Functions 


The test function, v, is assumed as in reference 14 
as 


vW = ju\ w) x\ W \x) + ju\ e) x\ e) (x) . 


(82) 
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In this paper, the test function components, Xi , are 
chosen as in the conventional MLPG method. The 

\4 

41 

►! 

« — l 

Z; W \x) components of the test functions are chosen as 

1 2 

9 

16 17 


power weight functions 14 , 


x\ W) (x) = { 



f dA 

2 " 

H 




\ S o J 



-i 4 


if 0 < </,■ < s 0 


(83) 


if dj > s 0 , 


with di = \\x - x,-| | . In Equation (83), s 0 is a user-defined 
parameter that determines the extent of the test 
functions (and hence Q s - see Figure 4). The 
components of the test functions chosen for 9 are the 
first derivatives of the components of the test functions 
chosen for the primary variable, w, i.e., 


(0) _ d Z, 


(w) 


Z = 


dx 


(84) 


as 9= ( dw/dx ) is also a primary variable. 


For this power function, the values of x\ , 
(dz\"^ 1 dx ) , and {d%\^ / dx) are zero when 
dj = S Q . As discussed in reference 14, when this test 
function is used, the k (nodc) in Equation (80d) reduces to 


. (node) 

K ij 


= El 


J 


-,(0 


dx 1 


dx~ 


•2 < e ) j 2 ,,,(»’) 

I d 4 dx 


-,(0 


dx z 


dx A 


\d^ytdx 


n(0 


dx 1 


dx 1 


j2 ( 0 ) J 2 I( .(0) 

| dz[ )d rj & 


n(0 


dx A 


dx z 


(85) 


Beam Configurations and Models 

A beam of constant flexural rigidity El and a length 
of 4/ is considered. The length 41 was specifically 
chosen to avoid scaling by a unit length, 1. Five models 
with 5, 9, 17, 33, and 65 nodes uniformly distributed 
along the length of the beam are considered. Figure 5 
shows a typical 17-node model. The distances between 
the nodes (Ax / I) in these models are 1, 0.5, 0.25, 
0.125, and 0.0625 for the 5-, 9-, 17-, 33-, and 65-node 
models, respectively. Numerical integration is used to 
integrate the system of equations as closed-form 
integration of the terms in Equations (80d) and (80f) is 
extremely complicated. 


Figure 5: A 17-node model of the beam 

Numerical Evaluations 


The radial basis MLPG (RPG) method was 
evaluated by applying the method to simple patch-test 
problems. The problems considered were (a) rigid 
body translation: 


w(x) = j3 0 , # = — = 0 , (86a) 

dx 

(b) rigid body rotation: 

w(x) = P\X, 9 = P i , (86b) 

and (c) constant-curvature condition: 

w(x) = /Gx 2 /2, 0 = P 2 X , (86c) 

where p 0 , P 1 , and /? 2 are arbitrary constants. The third 
patch test is equivalent to the problem of a cantilever 
beam with a moment, M=EI(d 2 w/dx 2 )= £7/%, applied at 
x=4/. The deflection, w, and the slope, 9, 
corresponding to problems (a), (b), and (c) were 
prescribed as essential boundary conditions (EBCs) at 
x=0 and x=4/. With these EBCs, the beam problems 
were analyzed using the RPG method with no 
polynomial basis. If the RPG method recovers the 
exact solution at all the interior nodes and at every 
arbitrary point of the beam, then the method passes the 
patch test. Note that in this work, near recovery of the 
exact solution is sufficient to pass the patch tests as the 
radial basis functions alone cannot represent 

polynomial solutions exactly. 

Compact RBF 

The compact radial functions described by 
Equations (36) and (37) were considered first. When 
using the compactly supported functions, the Equations 
(59) - (66) are evaluated with N = «, the number of 
nodes in the influence domain of the point x under 
consideration. 18 ' 19 This use of the compact functions 
forces the [Q B ] of Equation (59) to become a (2 n, 2 n) 
matrix that must be evaluated once for every node in 
the model. 
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The RPG method with no polynomial terms was 
unable to reproduce the exact solutions in Equation 
(86), and thus failed the patch tests. A quadratic 
polynomial basis (m = 3; p w : (1, x, x 2 ), pg. (0, 1, 2x)) 
was used in Equation (67), increasing the size of the 
[Q b ] matrix to (2n+in, 2n+m). The RPG method with 
polynomial basis (the hybrid RPG method) reproduced 
the solutions in Equation (86) to machine accuracy, 
thus passing the patch tests. 

Next, mixed boundary value problems were 
considered. The first problem considered was a 
cantilever beam with a tip load (Figure 6). Because the 
exact solution for this problem is cubic in x, the hybrid 
RPG method with a cubic polynomial basis function 
reproduced the exact solution. A simply supported 
beam subjected to a uniformly distributed load (Figure 
7) was considered next. Because the exact solution for 
this problem is quartic in x, the hybrid RPG method 
with quartic basis yielded the solution exactly. 

z, w 


t 



Figure 6: Cantilever beam with a tip load 


z, w 

t 



Figure 7: Simply supported beam subjected to a 
uniformly distributed load 


As with the finite element method and the 
conventional MLPG method 13, 14 , the hybrid RPG 
algorithm should be robust enough to yield good 
solutions when a low order polynomial basis function is 
used. A convergence test was conducted to study the 
performance of the hybrid method in solving the 
problems in Figures 6 and 7. A quadratic polynomial 
basis function was used. For all models (5, 9, 17, 33, 
and 65 nodes), the method did not yield meaningful 
results. Thus, it was concluded that as long as the order 
of the polynomial basis was sufficient to reproduce the 
solution exactly, the polynomial terms overpowered the 
radial basis functions. This condition is too restrictive, 
and hence compact radial functions are dropped from 
further consideration. 


Non-compactly Supported RBF 

Because the compactly supported radial basis 
functions are incapable of producing meaningful results 
for beam problems, the non-compact functions of Table 
1 are considered. In these functions, 

d i 

r = ^~, (87) 

s j 

where dj = ||x-x,-||, and Sj is some normalizing distance, 
usually chosen to be the entire problem domain, Q (in 
this work, 0 < x < L ). As ,v ; covers the entire problem 
domain, [Q B ] is an ( N , N) matrix that is evaluated and 
inverted once. 

Upon implementation of the functions in Table 1, 
the cubic function, 

R(r) = r\ (88) 

worked very well for the current C 1 problems. The 
RPG method with no polynomial basis and using 
Equation (88) was applied to the patch tests represented 
by Equations (86). The method successfully 
reproduced the exact solutions to machine accuracy, 
thus passing all the patch tests. Additionally, all 
functions of the form 

R(r) = r (2z_1) , (89) 

where z > 1, performed successfully, though r 3 gave the 
best results. 

Next, the RPG method with the RBF in Equation 
(88) was used to solve mixed boundary value problems. 
In the method, a 12-point Gaussian integration was 
used, the value of (s„ / I), which defines the extent of the 
test functions (see Equation (83)), was set as 
(s 0 / / =2 Ax), and the value of (sj / /), which defines the 
extent of the trial functions (Equation (87)), was set as 
(sj / l = L). For the cantilever beam with a tip load in 
Figure 6, the RPG method yielded excellent results. 
The simply supported beam problem with a uniformly 
distributed load (Figure 7) was analyzed using 17-, 33-, 
and 65-node models. The maximum deflection values, 
i.e., the deflection at (x = L / 2), for these three models 
obtained using the RPG method and using the 
conventional MLPG method with a quadratic 
polynomial basis function are compared in Table 2. In 
the MLPG method, a 20-point Gaussian integration was 
used, the value of (.? 0 / /) was set as (.? 0 / / =2Ax), and the 
value of (sj / l) was set as (sj / / =8 At). From this table, 
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it is seen that the RPG method performs as accurately 
as the conventional MLPG method. 


the perceived advantages of the RPG method over the 
MLPG method. 


Table 2: Maximum deflection values for three nodal 
models obtained using the RPG and MLPG methods 
compared to the exact solution. 


Model 

Maximum deflection (at x = L / 2) 

Exact 

RPG 

MLPG 

17-node 

-3.3333e-7 

-3.2739e-7 

-3.3106e-7 

3 3 -node 

-3.3333e-7 

-3.3407e-7 

-3.3735e-7 

65-node 

-3.3333e-7 

-3.3420e-7 

-3.3848e-7 




Exact solutions 


O MtPG 
X A^mlpg 

□ Frpg 


(b) Secondary variables (NOTE: P MLPG not shown) 

Figure 8: RPG, MLPG, and Exact solutions obtained 
using the 65-node model for the simply supported beam 
subjected to a uniformly distributed load 

The results obtained for deflection, slope, moment, 
and shear using the 65-node model are presented in 
Figure 8. In this figure, the RPG results are compared 
to the exact solution and to the solution obtained using 
the conventional MLPG method with a quadratic 
polynomial basis function. For each of the nodal 
models (17, 33, and 65 nodes), the RPG values for 
deflection, slope, and moment were as accurate as the 
MLPG values and were in excellent agreement with the 
exact values. In addition, the RPG values for shear 
converged with model refinement. The MLPG solution 
for the shear was erratic, and is not shown in Figure 8. 
The quadratic basis function is insufficient to accurately 
calculate the third derivatives for this problem, and the 
method could not recover the values with model 
refinement; the solution for the shear converged only as 
the order of the basis function was increased to 
quartic. 13 The results discussed for this problem verify 


The RPG method with the RBF in Equation (88) 
was then applied to a problem with load discontinuity. 
The problem considered was the cantilever beam with a 
uniformly distributed load on a portion of the beam 
shown in Figure 9. The RPG solution (with 
(s 0 / 1= 4Ax)) for the cantilever beam problem exhibited 
convergence with model refinement. These results are 
consistent with those reported in reference 14, where 
this problem was studied using the conventional MLPG 
method. The exact, MLPG, and RPG values for 
deflection and moment for this problem obtained using 
a 65-node model are compared in Figure 10. The 
parameters used for the MLPG method are the same as 
those reported above for the simply supported beam 
problem. The RPG method handled the load 
discontinuity well and yielded results in overall 
agreement with the exact solutions. 


/* 


i r 


i 


Figure 9: Cantilever beam with a uniformly distributed 
load on a portion of the beam 



Exact solutions 


Figure 10: RPG, MLPG, and Exact solutions obtained 
using the 65-node model for the cantilever beam with a 
uniformly distributed load on a portion of the beam 


Concluding Remarks 

A radial basis function implementation of the 
MLPG method was presented to study Euler-Bernoulli 
beam problems. Like the conventional MLPG method, 
this radial basis variation (RPG) is based on the local 
weak form developed from the classical weighted 
residual form of the fourth-order governing differential 
equation. In this method, radial basis functions, rather 
than generalized moving least squares (GMLS) 
interpolations, were used to develop the trial functions, 
and test functions were chosen as simple weight 
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functions as in the conventional MLPG method. RPG 
equations were developed with and without including 
polynomial basis function terms. 

The compactly supported radial basis functions did 
not perform well without polynomial terms in the 
computations. When polynomial terms were included, 
the compactly supported RPG method passed the patch 
tests. However, the method did not yield meaningful 
results for mixed boundary value problems unless the 
order of the polynomial basis function was of the same 
order as the exact solution of the problem. This result 
restricts the use of the method. The use of compactly 
supported radial basis functions is not recommended. 

The non-compactly supported cubic radial basis 
function performed very well when no polynomial 
terms were included in the computations. The RPG 
method with the cubic radial basis function passed all 
the patch tests and yielded results for mixed boundary 
value problems that are comparable to those obtained 
using the conventional MLPG method. The RPG 
method with a cubic radial basis also yielded very good 
results for a problem with discontinuous loading 
conditions. The accuracy of solutions obtained by the 
RPG method, combined with the computational 
efficiency of using the radial basis functions rather than 
the GMLS interpolations to approximate the trial 
functions, makes the RPG method a very attractive 
variation of the MLPG method. 
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